library(tidyverse)
library(here)
library(vegan)
library(RColorBrewer)
library(broom)
source(here("util_functions.R"))

Read data

criteria_path <- here("data/disease_criteria")
parse_diagnoses <- function(sim_results){
  
  tibble(criteria = rownames(sim_results)) %>% 
    mutate(data = map(criteria, function(c){
      do.call(rbind, sim_results[c,]) %>% 
        mutate(diagnosis = map(diagnoses, process_diagnoses)) %>% 
        unnest(diagnosis) %>% 
        select(diagnosis)
      })) %>% 
    unnest(data)
         
}
# Read in ChatGPT output data for additional conditions and parse
df <- tibble(file = list.files(here("data/chatgpt_output/other_conditions"), full.names = T))%>% 
  mutate(data = map(file, readRDS)) %>% 
  mutate(data = map(data, parse_diagnoses)) %>% 
  unnest(data) %>% 
  select(-file)

# Read in original MCAS chatGPT output and append to data from additional conditions
df <- read_csv(here("data/compiled_mcas_chatgpt_diagnoses.csv")) %>% 
  mutate(criteria = gsub("consensus", "mcas_", consensus)) %>% 
  select(criteria, diagnosis) %>% 
  bind_rows(df) 
Rows: 38029 Columns: 2── Column specification ─────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): diagnosis, consensus
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
df
# write_csv(df, here("data/compiled_all_chatgpt_diagnoses.csv"))
df <- read_csv(here("data/compiled_all_chatgpt_diagnoses.csv"))
Rows: 117354 Columns: 2── Column specification ─────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): criteria, diagnosis
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
b <- df %>% count(diagnosis)

Plot distributions

# Rank abundance plot
custom_pal <- RColorBrewer::brewer.pal(7, "Set1")[-6]

df %>% 
  count(criteria, diagnosis, sort = T) %>% 
  group_by(criteria) %>% 
  mutate(freq = n/sum(n)) %>% 
  arrange(desc(freq), .by_group = T) %>% 
  mutate(rank = 1:n()) %>% 
  select(criteria, freq, rank) %>% 
  mutate(cs = cumsum(freq)) %>% 
  filter(rank <= 50) %>% 
  mutate(criteria = toupper(criteria)) %>% 
  ggplot(aes(x = rank, y = freq, color = factor(criteria)))+
    geom_line() +
    theme_classic()+
    scale_color_manual(values = custom_pal)+
  scale_x_continuous(expand = c(0,0))+
  scale_y_continuous(expand = c(0,0)) +
  labs(x = "Diagnosis rank", y = "Diagnosis frequency", color = "")

Calculate and plot diversity

# Calculate diversity values for diagnosis distributions
diagnosis_table <- table(df$criteria, df$diagnosis)
div_diag  <- vegan::diversity(diagnosis_table)
div_diag  
 aha_kawasaki eular_acr_sle        mcas_1        mcas_2      migraine     slicc_sle 
     4.768183      4.749758      4.847143      6.295376      4.134877      4.393290 
# Calculate similarity
plot_similarity <- function(d_method){
  broom::tidy(1-vegan::vegdist(diagnosis_table, method = d_method)) %>% 
    mutate(item = sprintf("%s_%s", item1, item2)) %>% 
    ggplot(aes(x = item, y = distance)) +
    geom_bar(stat = "identity")+
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
    coord_flip() +
    labs(x="", y="Similarity")
}

plot_similarity("bray")

(1-as.matrix(vegan::vegdist(diagnosis_table, method = "bray"))) %>%
  ComplexHeatmap::Heatmap(.,
                          name = "Similarity",
                          col = viridis::viridis(100))

diagnosis_pca <- as.data.frame(diagnosis_table) %>% 
  pivot_wider(names_from = "Var2", values_from = "Freq", values_fill = 0) %>% 
  column_to_rownames("Var1") %>% 
  prcomp()

as.data.frame(diagnosis_pca$x) %>% 
  rownames_to_column("criteria") %>% 
  ggplot(aes(x = PC1, y = PC2, label = criteria))+
  geom_label() +
  theme_bw()

# Plot diversity values
broom::tidy(div_diag) %>% 
  mutate(names = toupper(names)) %>% 
  ggplot(aes(x=names, y = x))+
  geom_bar(stat = "identity") +
  theme_bw() +
  labs(x = "", y = "Shannon diversity") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
Warning: 'tidy.numeric' is deprecated.
See help("Deprecated")

Statistical testing

All group randomization

pairwise_div_difference <- function(df){
  df_div <- enframe(vegan::diversity(table(df$criteria, df$diagnosis)))
  df_diff <- data.frame(t(combn(unique(sort(df$criteria)),2))) %>% 
    left_join(df_div, by = c("X1"="name")) %>% 
    left_join(df_div, by = c("X2"="name")) %>% 
    unite("pair", X1, X2, sep = ".") %>%
    mutate(entropy_difference = value.x - value.y) %>% 
    select(-contains("value"))
  deframe(df_diff)
}

# pairwise_div_difference(df)
enframe(pairwise_div_difference(df), name = "pair", value = "entropy_difference")
pairwise_div_iteration <- function(){
  df %>% 
    count(criteria) %>% 
    mutate(diagnosis = map(n, function(x) sample(df$diagnosis, x, replace = T))) %>% 
    select(-n) %>% 
    unnest(diagnosis) %>% 
    pairwise_div_difference(.)
}

pairwise_div_iteration()
aha_kawasaki.eular_acr_sle        aha_kawasaki.mcas_1        aha_kawasaki.mcas_2 
              -0.008694848                0.008252266                0.012165368 
     aha_kawasaki.migraine     aha_kawasaki.slicc_sle       eular_acr_sle.mcas_1 
               0.002825259               -0.001106212                0.016947114 
      eular_acr_sle.mcas_2     eular_acr_sle.migraine    eular_acr_sle.slicc_sle 
               0.020860217                0.011520107                0.007588636 
             mcas_1.mcas_2            mcas_1.migraine           mcas_1.slicc_sle 
               0.003913103               -0.005427007               -0.009358478 
           mcas_2.migraine           mcas_2.slicc_sle         migraine.slicc_sle 
              -0.009340110               -0.013271581               -0.003931471 
set.seed(1234)
pairs_permutation_distributions <- data.frame(t(replicate(10000, pairwise_div_iteration())))
[WARNING] Deprecated: --self-contained. use --embed-resources --standalone
sapply(combn(unique(sort(df$criteria)),2, simplify = F), paste, collapse = ".")
 [1] "aha_kawasaki.eular_acr_sle" "aha_kawasaki.mcas_1"        "aha_kawasaki.mcas_2"       
 [4] "aha_kawasaki.migraine"      "aha_kawasaki.slicc_sle"     "eular_acr_sle.mcas_1"      
 [7] "eular_acr_sle.mcas_2"       "eular_acr_sle.migraine"     "eular_acr_sle.slicc_sle"   
[10] "mcas_1.mcas_2"              "mcas_1.migraine"            "mcas_1.slicc_sle"          
[13] "mcas_2.migraine"            "mcas_2.slicc_sle"           "migraine.slicc_sle"        
names(pairs_permutation_distributions)
 [1] "aha_kawasaki.eular_acr_sle" "aha_kawasaki.mcas_1"        "aha_kawasaki.mcas_2"       
 [4] "aha_kawasaki.migraine"      "aha_kawasaki.slicc_sle"     "eular_acr_sle.mcas_1"      
 [7] "eular_acr_sle.mcas_2"       "eular_acr_sle.migraine"     "eular_acr_sle.slicc_sle"   
[10] "mcas_1.mcas_2"              "mcas_1.migraine"            "mcas_1.slicc_sle"          
[13] "mcas_2.migraine"            "mcas_2.slicc_sle"           "migraine.slicc_sle"        
pairs_entropy_diff <- enframe(pairwise_div_difference(df), name = "pair", value = "entropy_difference") %>% 
  mutate(entropy_difference = ifelse(entropy_difference >=0, 0-entropy_difference, entropy_difference))
pairs_entropy_diff
read_csv(here("data/all_criteria_null_difference_distributions.csv")) %>% 
  # pivot_longer(everything(), names_to = "pair", values_to = "entropy_difference") %>% 
  ggplot(aes(x = entropy_difference)) +
  geom_histogram(binwidth = 0.01)+
  geom_vline(data = pairs_entropy_diff, aes(xintercept = entropy_difference), color = "red")+
  facet_wrap(~pair, scales = "free")+
  theme_classic()
Rows: 150000 Columns: 2── Column specification ─────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): pair
dbl (1): entropy_difference
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Calculate p-values for diversity differences based on the quantile of the 
# diversity difference value compared to the null hypothesis diversity distribution values
# Pairwise comparison adjusted
entropy_diff_table <- read_csv(here("data/all_criteria_null_difference_distributions.csv"))  %>% 
  # pivot_longer(everything(), names_to = "pair", values_to = "entropy_difference") %>% 
  group_by(pair) %>% 
  summarise(entropy_difference_distribution = list(entropy_difference)) %>% 
  right_join(pairs_entropy_diff) %>% 
  mutate(p_value = map2_dbl(entropy_difference, entropy_difference_distribution, ~ecdf(.y)(.x))) %>% 
  mutate(p_adj = p.adjust(p_value, method = "bonferroni", n = n()))
Rows: 150000 Columns: 2── Column specification ─────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): pair
dbl (1): entropy_difference
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Joining with `by = join_by(pair)`
entropy_diff_table
# Plot diversity values
broom::tidy(div_diag) %>% 
  mutate(names = tolower(names)) %>% 
  ggplot(aes(x=names, y = x))+
  geom_bar(stat = "identity") +
  ggpubr::stat_pvalue_manual(
    entropy_diff_table %>% separate(pair, into = c("group1", "group2"), sep = "\\."),
    label= "p_adj",
    y.position = 7, step.increase = 0.2
  ) +
  theme_bw() +
  labs(x = "", y = "Shannon diversity") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
Warning: 'tidy.numeric' is deprecated.
See help("Deprecated")

Pairwise randomization

pairs_list <- as.list(as.data.frame(combn(unique(df$criteria),2)))
names(pairs_list) <- sapply(pairs_list, paste, collapse = ".")
pairs_list
$mcas_1.mcas_2
[1] "mcas_1" "mcas_2"

$mcas_1.aha_kawasaki
[1] "mcas_1"       "aha_kawasaki"

$mcas_1.eular_acr_sle
[1] "mcas_1"        "eular_acr_sle"

$mcas_1.migraine
[1] "mcas_1"   "migraine"

$mcas_1.slicc_sle
[1] "mcas_1"    "slicc_sle"

$mcas_2.aha_kawasaki
[1] "mcas_2"       "aha_kawasaki"

$mcas_2.eular_acr_sle
[1] "mcas_2"        "eular_acr_sle"

$mcas_2.migraine
[1] "mcas_2"   "migraine"

$mcas_2.slicc_sle
[1] "mcas_2"    "slicc_sle"

$aha_kawasaki.eular_acr_sle
[1] "aha_kawasaki"  "eular_acr_sle"

$aha_kawasaki.migraine
[1] "aha_kawasaki" "migraine"    

$aha_kawasaki.slicc_sle
[1] "aha_kawasaki" "slicc_sle"   

$eular_acr_sle.migraine
[1] "eular_acr_sle" "migraine"     

$eular_acr_sle.slicc_sle
[1] "eular_acr_sle" "slicc_sle"    

$migraine.slicc_sle
[1] "migraine"  "slicc_sle"
# Create function that subsets df of all diagnoses for a criteria based on a pair of two criteria to compare
# Then randomizes the diagnoses across the two criteria, 
# measures the resulting diversity of each criteria's diagnoses
# and calculates a difference in the diversity value
# randomize can be set F to obtain original difference in diversity w/o randomization
# To be used with replicate to generate a null hypothesis distribution 
permutation_pair_iteration <- function(pair, randomize = T){
  df <- df %>% 
    filter(criteria %in% pair) 
  if (randomize == T){
    df <- df %>% mutate(diagnosis = sample(diagnosis, n()))
  }
  div <- vegan::diversity(table(df$criteria, df$diagnosis))
  unname(div[1] - div[2])
}

replicate(10,permutation_pair_iteration(pairs_list[[1]], randomize = T))
 [1]  0.0049722988  0.0081672943  0.0156799259  0.0038521584  0.0060045361 -0.0151564640
 [7]  0.0144935438 -0.0333420796  0.0278905935 -0.0006473075
replicate(10,permutation_pair_iteration(pairs_list[[1]], randomize = F))
 [1] -1.434911 -1.434911 -1.434911 -1.434911 -1.434911 -1.434911 -1.434911 -1.434911 -1.434911
[10] -1.434911
# # Create a null hypothesis distribution for the difference in diversity
# # For all pairs of criteria

# set.seed(1234)
# pairs_permutation_distributions <- lapply(pairs_list, function(x) replicate(10000, permutation_pair_iteration(x)))
# names(pairs_permutation_distributions) <- lapply(pairs_list, function(x) paste(x, collapse = "-"))
# data.frame(pairs_permutation_distributions)
# write_csv(data.frame(pairs_permutation_distributions), here("data/criteria_pairwise_null_difference_distributions.csv"))

pairs_permutation_distributions <- read_csv(here("data/criteria_pairwise_null_difference_distributions.csv"))
Rows: 10000 Columns: 15── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (15): mcas_1.mcas_2, mcas_1.aha_kawasaki, mcas_1.eular_acr_sle, mcas_1.migraine, mcas_1.slicc_sle, mcas_2.aha_kawasaki, mcas_2.eular_acr_sle, mcas_2.migraine, mcas_2.slicc_sle, aha_kawasaki.eular_...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pairs_permutation_distributions
# Create a df of the actual difference in diversity between pairs of criteria
pairs_entropy_diff <- lapply(pairs_list, permutation_pair_iteration, randomize = F)
names(pairs_entropy_diff) <- lapply(pairs_list, function(x) paste(x, collapse = "-"))
pairs_entropy_diff <- data.frame(pairs_entropy_diff) %>%  
  pivot_longer(everything(), names_to = "pair", values_to = "entropy_difference") %>% 
  mutate(entropy_difference = ifelse(entropy_difference >=0, 0-entropy_difference, entropy_difference))
pairs_entropy_diff
# Create plots of null hypothesis distributions, overlaid with actual diversity difference
read_csv(here("data/all_criteria_pairwise_null_difference_distributions.csv")) %>%
  group_by(pair) %>% 
  nest() %>% 
  separate(pair, into = c("p1","p2"), sep = "\\.", remove = F) %>%
  rowwise() %>% mutate(pair=paste(sort(c(p1,p2)), collapse = ".")) %>% 
  select(!pair:p2) %>% 
  unnest(data) %>% 
  # pivot_longer(cols = everything(), names_to = "pair", values_to = "entropy_difference") %>% 
  ggplot(aes(x = entropy_difference)) +
  geom_histogram(binwidth = 0.01)+
  geom_vline(data = pairs_entropy_diff, aes(xintercept = entropy_difference), color = "red")+
  facet_wrap(~pair, scales = "free")+
  theme_classic()
Rows: 150000 Columns: 2── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): pair
dbl (1): entropy_difference
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Adding missing grouping variables: `pair`

# Calculate p-values for diversity differences based on the quantile of the 
# diversity difference value compared to the null hypothesis diversity distribution values
# Pairwise comparison adjusted
data.frame(pairs_permutation_distributions) %>% 
  pivot_longer(everything(), names_to = "pair", values_to = "entropy_difference") %>% 
  group_by(pair) %>% 
  summarise(entropy_difference_distribution = list(entropy_difference)) %>% 
  right_join(pairs_entropy_diff) %>% 
  mutate(p_value = map2_dbl(entropy_difference, entropy_difference_distribution, ~ecdf(.y)(.x))) %>% 
  mutate(p_adj = p.adjust(p_value, method = "bonferroni", n = n()))
Joining with `by = join_by(pair)`

Plots

LS0tCnRpdGxlOiAiQ29tcGFyaW5nIE1DQVMgY3JpdGVyaWEgdG8gb3RoZXIgY29uZGl0aW9ucyIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KGhlcmUpCmxpYnJhcnkodmVnYW4pCmxpYnJhcnkoUkNvbG9yQnJld2VyKQpsaWJyYXJ5KGJyb29tKQpzb3VyY2UoaGVyZSgidXRpbF9mdW5jdGlvbnMuUiIpKQpgYGAKCiMjIFJlYWQgZGF0YQoKYGBge3J9CmNyaXRlcmlhX3BhdGggPC0gaGVyZSgiZGF0YS9kaXNlYXNlX2NyaXRlcmlhIikKYGBgCgpgYGB7cn0KcGFyc2VfZGlhZ25vc2VzIDwtIGZ1bmN0aW9uKHNpbV9yZXN1bHRzKXsKICAKICB0aWJibGUoY3JpdGVyaWEgPSByb3duYW1lcyhzaW1fcmVzdWx0cykpICU+JSAKICAgIG11dGF0ZShkYXRhID0gbWFwKGNyaXRlcmlhLCBmdW5jdGlvbihjKXsKICAgICAgZG8uY2FsbChyYmluZCwgc2ltX3Jlc3VsdHNbYyxdKSAlPiUgCiAgICAgICAgbXV0YXRlKGRpYWdub3NpcyA9IG1hcChkaWFnbm9zZXMsIHByb2Nlc3NfZGlhZ25vc2VzKSkgJT4lIAogICAgICAgIHVubmVzdChkaWFnbm9zaXMpICU+JSAKICAgICAgICBzZWxlY3QoZGlhZ25vc2lzKQogICAgICB9KSkgJT4lIAogICAgdW5uZXN0KGRhdGEpCiAgICAgICAgIAp9CmBgYAoKYGBge3J9CiMgUmVhZCBpbiBDaGF0R1BUIG91dHB1dCBkYXRhIGZvciBhZGRpdGlvbmFsIGNvbmRpdGlvbnMgYW5kIHBhcnNlCmRmIDwtIHRpYmJsZShmaWxlID0gbGlzdC5maWxlcyhoZXJlKCJkYXRhL2NoYXRncHRfb3V0cHV0L290aGVyX2NvbmRpdGlvbnMiKSwgZnVsbC5uYW1lcyA9IFQpKSU+JSAKICBtdXRhdGUoZGF0YSA9IG1hcChmaWxlLCByZWFkUkRTKSkgJT4lIAogIG11dGF0ZShkYXRhID0gbWFwKGRhdGEsIHBhcnNlX2RpYWdub3NlcykpICU+JSAKICB1bm5lc3QoZGF0YSkgJT4lIAogIHNlbGVjdCgtZmlsZSkKCiMgUmVhZCBpbiBvcmlnaW5hbCBNQ0FTIGNoYXRHUFQgb3V0cHV0IGFuZCBhcHBlbmQgdG8gZGF0YSBmcm9tIGFkZGl0aW9uYWwgY29uZGl0aW9ucwpkZiA8LSByZWFkX2NzdihoZXJlKCJkYXRhL2NvbXBpbGVkX21jYXNfY2hhdGdwdF9kaWFnbm9zZXMuY3N2IikpICU+JSAKICBtdXRhdGUoY3JpdGVyaWEgPSBnc3ViKCJjb25zZW5zdXMiLCAibWNhc18iLCBjb25zZW5zdXMpKSAlPiUgCiAgc2VsZWN0KGNyaXRlcmlhLCBkaWFnbm9zaXMpICU+JSAKICBiaW5kX3Jvd3MoZGYpIAoKZGYKYGBgCgpgYGB7cn0KIyB3cml0ZV9jc3YoZGYsIGhlcmUoImRhdGEvY29tcGlsZWRfYWxsX2NoYXRncHRfZGlhZ25vc2VzLmNzdiIpKQpkZiA8LSByZWFkX2NzdihoZXJlKCJkYXRhL2NvbXBpbGVkX2FsbF9jaGF0Z3B0X2RpYWdub3Nlcy5jc3YiKSkKYGBgCgoKIyMgUGxvdCBkaXN0cmlidXRpb25zIAoKYGBge3J9CiMgUmFuayBhYnVuZGFuY2UgcGxvdApjdXN0b21fcGFsIDwtIFJDb2xvckJyZXdlcjo6YnJld2VyLnBhbCg3LCAiU2V0MSIpWy02XQoKZGYgJT4lIAogIGNvdW50KGNyaXRlcmlhLCBkaWFnbm9zaXMsIHNvcnQgPSBUKSAlPiUgCiAgZ3JvdXBfYnkoY3JpdGVyaWEpICU+JSAKICBtdXRhdGUoZnJlcSA9IG4vc3VtKG4pKSAlPiUgCiAgYXJyYW5nZShkZXNjKGZyZXEpLCAuYnlfZ3JvdXAgPSBUKSAlPiUgCiAgbXV0YXRlKHJhbmsgPSAxOm4oKSkgJT4lIAogIHNlbGVjdChjcml0ZXJpYSwgZnJlcSwgcmFuaykgJT4lIAogIG11dGF0ZShjcyA9IGN1bXN1bShmcmVxKSkgJT4lIAogIGZpbHRlcihyYW5rIDw9IDUwKSAlPiUgCiAgbXV0YXRlKGNyaXRlcmlhID0gdG91cHBlcihjcml0ZXJpYSkpICU+JSAKICBnZ3Bsb3QoYWVzKHggPSByYW5rLCB5ID0gZnJlcSwgY29sb3IgPSBmYWN0b3IoY3JpdGVyaWEpKSkrCiAgICBnZW9tX2xpbmUoKSArCiAgICB0aGVtZV9jbGFzc2ljKCkrCiAgICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzID0gY3VzdG9tX3BhbCkrCiAgc2NhbGVfeF9jb250aW51b3VzKGV4cGFuZCA9IGMoMCwwKSkrCiAgc2NhbGVfeV9jb250aW51b3VzKGV4cGFuZCA9IGMoMCwwKSkgKwogIGxhYnMoeCA9ICJEaWFnbm9zaXMgcmFuayIsIHkgPSAiRGlhZ25vc2lzIGZyZXF1ZW5jeSIsIGNvbG9yID0gIiIpCmBgYAojIyBDYWxjdWxhdGUgYW5kIHBsb3QgZGl2ZXJzaXR5CgpgYGB7cn0KIyBDYWxjdWxhdGUgZGl2ZXJzaXR5IHZhbHVlcyBmb3IgZGlhZ25vc2lzIGRpc3RyaWJ1dGlvbnMKZGlhZ25vc2lzX3RhYmxlIDwtIHRhYmxlKGRmJGNyaXRlcmlhLCBkZiRkaWFnbm9zaXMpCmRpdl9kaWFnICA8LSB2ZWdhbjo6ZGl2ZXJzaXR5KGRpYWdub3Npc190YWJsZSkKZGl2X2RpYWcgIApgYGAKCmBgYHtyfQojIENhbGN1bGF0ZSBzaW1pbGFyaXR5CnBsb3Rfc2ltaWxhcml0eSA8LSBmdW5jdGlvbihkX21ldGhvZCl7CiAgYnJvb206OnRpZHkoMS12ZWdhbjo6dmVnZGlzdChkaWFnbm9zaXNfdGFibGUsIG1ldGhvZCA9IGRfbWV0aG9kKSkgJT4lIAogICAgbXV0YXRlKGl0ZW0gPSBzcHJpbnRmKCIlc18lcyIsIGl0ZW0xLCBpdGVtMikpICU+JSAKICAgIGdncGxvdChhZXMoeCA9IGl0ZW0sIHkgPSBkaXN0YW5jZSkpICsKICAgIGdlb21fYmFyKHN0YXQgPSAiaWRlbnRpdHkiKSsKICAgIHRoZW1lX2J3KCkgKwogICAgdGhlbWUoYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSA5MCwgaGp1c3QgPSAxKSkgKwogICAgY29vcmRfZmxpcCgpICsKICAgIGxhYnMoeD0iIiwgeT0iU2ltaWxhcml0eSIpCn0KCnBsb3Rfc2ltaWxhcml0eSgiYnJheSIpCmBgYApgYGB7cn0KIyBTaW1pbGFyaXR5IGhlYXRtYXAKKDEtYXMubWF0cml4KHZlZ2FuOjp2ZWdkaXN0KGRpYWdub3Npc190YWJsZSwgbWV0aG9kID0gImJyYXkiKSkpICU+JQogIENvbXBsZXhIZWF0bWFwOjpIZWF0bWFwKC4sCiAgICAgICAgICAgICAgICAgICAgICAgICAgbmFtZSA9ICJTaW1pbGFyaXR5IiwKICAgICAgICAgICAgICAgICAgICAgICAgICBjb2wgPSB2aXJpZGlzOjp2aXJpZGlzKDEwMCkpCmBgYAoKCmBgYHtyfQojIFNpbWlsYXJpdHkgUENBIApkaWFnbm9zaXNfcGNhIDwtIGFzLmRhdGEuZnJhbWUoZGlhZ25vc2lzX3RhYmxlKSAlPiUgCiAgcGl2b3Rfd2lkZXIobmFtZXNfZnJvbSA9ICJWYXIyIiwgdmFsdWVzX2Zyb20gPSAiRnJlcSIsIHZhbHVlc19maWxsID0gMCkgJT4lIAogIGNvbHVtbl90b19yb3duYW1lcygiVmFyMSIpICU+JSAKICBwcmNvbXAoKQoKYXMuZGF0YS5mcmFtZShkaWFnbm9zaXNfcGNhJHgpICU+JSAKICByb3duYW1lc190b19jb2x1bW4oImNyaXRlcmlhIikgJT4lIAogIGdncGxvdChhZXMoeCA9IFBDMSwgeSA9IFBDMiwgbGFiZWwgPSBjcml0ZXJpYSkpKwogIGdlb21fbGFiZWwoKSArCiAgdGhlbWVfYncoKQpgYGAKCmBgYHtyLCBmaWcud2lkdGggPSA1LCBmaWcuaGVpZ2h0PTV9CiMgUGxvdCBkaXZlcnNpdHkgdmFsdWVzCmJyb29tOjp0aWR5KGRpdl9kaWFnKSAlPiUgCiAgbXV0YXRlKG5hbWVzID0gdG91cHBlcihuYW1lcykpICU+JSAKICBnZ3Bsb3QoYWVzKHg9bmFtZXMsIHkgPSB4KSkrCiAgZ2VvbV9iYXIoc3RhdCA9ICJpZGVudGl0eSIpICsKICB0aGVtZV9idygpICsKICBsYWJzKHggPSAiIiwgeSA9ICJTaGFubm9uIGRpdmVyc2l0eSIpICsKICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZSA9IDQ1LCBoanVzdCA9IDEpKQpgYGAKCiMjIFN0YXRpc3RpY2FsIHRlc3RpbmcKCiMjIyBBbGwgZ3JvdXAgcmFuZG9taXphdGlvbgoKYGBge3J9CnBhaXJ3aXNlX2Rpdl9kaWZmZXJlbmNlIDwtIGZ1bmN0aW9uKGRmKXsKICBkZl9kaXYgPC0gZW5mcmFtZSh2ZWdhbjo6ZGl2ZXJzaXR5KHRhYmxlKGRmJGNyaXRlcmlhLCBkZiRkaWFnbm9zaXMpKSkKICBkZl9kaWZmIDwtIGRhdGEuZnJhbWUodChjb21ibih1bmlxdWUoc29ydChkZiRjcml0ZXJpYSkpLDIpKSkgJT4lIAogICAgbGVmdF9qb2luKGRmX2RpdiwgYnkgPSBjKCJYMSI9Im5hbWUiKSkgJT4lIAogICAgbGVmdF9qb2luKGRmX2RpdiwgYnkgPSBjKCJYMiI9Im5hbWUiKSkgJT4lIAogICAgdW5pdGUoInBhaXIiLCBYMSwgWDIsIHNlcCA9ICIuIikgJT4lCiAgICBtdXRhdGUoZW50cm9weV9kaWZmZXJlbmNlID0gdmFsdWUueCAtIHZhbHVlLnkpICU+JSAKICAgIHNlbGVjdCgtY29udGFpbnMoInZhbHVlIikpCiAgZGVmcmFtZShkZl9kaWZmKQp9CgojIHBhaXJ3aXNlX2Rpdl9kaWZmZXJlbmNlKGRmKQplbmZyYW1lKHBhaXJ3aXNlX2Rpdl9kaWZmZXJlbmNlKGRmKSwgbmFtZSA9ICJwYWlyIiwgdmFsdWUgPSAiZW50cm9weV9kaWZmZXJlbmNlIikKYGBgCgpgYGB7cn0KcGFpcndpc2VfZGl2X2l0ZXJhdGlvbiA8LSBmdW5jdGlvbigpewogIGRmICU+JSAKICAgIGNvdW50KGNyaXRlcmlhKSAlPiUgCiAgICBtdXRhdGUoZGlhZ25vc2lzID0gbWFwKG4sIGZ1bmN0aW9uKHgpIHNhbXBsZShkZiRkaWFnbm9zaXMsIHgsIHJlcGxhY2UgPSBUKSkpICU+JSAKICAgIHNlbGVjdCgtbikgJT4lIAogICAgdW5uZXN0KGRpYWdub3NpcykgJT4lIAogICAgcGFpcndpc2VfZGl2X2RpZmZlcmVuY2UoLikKfQoKcGFpcndpc2VfZGl2X2l0ZXJhdGlvbigpCmBgYAoKYGBge3J9CnNldC5zZWVkKDEyMzQpCnBhaXJzX3Blcm11dGF0aW9uX2Rpc3RyaWJ1dGlvbnMgPC0gZGF0YS5mcmFtZSh0KHJlcGxpY2F0ZSgxMDAwMCwgcGFpcndpc2VfZGl2X2l0ZXJhdGlvbigpKSkpCmBgYApgYGB7cn0Kc2FwcGx5KGNvbWJuKHVuaXF1ZShzb3J0KGRmJGNyaXRlcmlhKSksMiwgc2ltcGxpZnkgPSBGKSwgcGFzdGUsIGNvbGxhcHNlID0gIi4iKQpuYW1lcyhwYWlyc19wZXJtdXRhdGlvbl9kaXN0cmlidXRpb25zKQpgYGAKYGBge3J9CnBhaXJzX2VudHJvcHlfZGlmZiA8LSBlbmZyYW1lKHBhaXJ3aXNlX2Rpdl9kaWZmZXJlbmNlKGRmKSwgbmFtZSA9ICJwYWlyIiwgdmFsdWUgPSAiZW50cm9weV9kaWZmZXJlbmNlIikgJT4lIAogIG11dGF0ZShlbnRyb3B5X2RpZmZlcmVuY2UgPSBpZmVsc2UoZW50cm9weV9kaWZmZXJlbmNlID49MCwgMC1lbnRyb3B5X2RpZmZlcmVuY2UsIGVudHJvcHlfZGlmZmVyZW5jZSkpCnBhaXJzX2VudHJvcHlfZGlmZgpgYGAKCgpgYGB7ciwgZmlnLndpZHRoPTEyLCBmaWcuaGVpZ2h0ID0gOH0KcmVhZF9jc3YoaGVyZSgiZGF0YS9hbGxfY3JpdGVyaWFfbnVsbF9kaWZmZXJlbmNlX2Rpc3RyaWJ1dGlvbnMuY3N2IikpICU+JSAKICAjIHBpdm90X2xvbmdlcihldmVyeXRoaW5nKCksIG5hbWVzX3RvID0gInBhaXIiLCB2YWx1ZXNfdG8gPSAiZW50cm9weV9kaWZmZXJlbmNlIikgJT4lIAogIGdncGxvdChhZXMoeCA9IGVudHJvcHlfZGlmZmVyZW5jZSkpICsKICBnZW9tX2hpc3RvZ3JhbShiaW53aWR0aCA9IDAuMDEpKwogIGdlb21fdmxpbmUoZGF0YSA9IHBhaXJzX2VudHJvcHlfZGlmZiwgYWVzKHhpbnRlcmNlcHQgPSBlbnRyb3B5X2RpZmZlcmVuY2UpLCBjb2xvciA9ICJyZWQiKSsKICBmYWNldF93cmFwKH5wYWlyLCBzY2FsZXMgPSAiZnJlZSIpKwogIHRoZW1lX2NsYXNzaWMoKQpgYGAKCmBgYHtyfQojIENhbGN1bGF0ZSBwLXZhbHVlcyBmb3IgZGl2ZXJzaXR5IGRpZmZlcmVuY2VzIGJhc2VkIG9uIHRoZSBxdWFudGlsZSBvZiB0aGUgCiMgZGl2ZXJzaXR5IGRpZmZlcmVuY2UgdmFsdWUgY29tcGFyZWQgdG8gdGhlIG51bGwgaHlwb3RoZXNpcyBkaXZlcnNpdHkgZGlzdHJpYnV0aW9uIHZhbHVlcwojIFBhaXJ3aXNlIGNvbXBhcmlzb24gYWRqdXN0ZWQKZW50cm9weV9kaWZmX3RhYmxlIDwtIHJlYWRfY3N2KGhlcmUoImRhdGEvYWxsX2NyaXRlcmlhX251bGxfZGlmZmVyZW5jZV9kaXN0cmlidXRpb25zLmNzdiIpKSAgJT4lIAogICMgcGl2b3RfbG9uZ2VyKGV2ZXJ5dGhpbmcoKSwgbmFtZXNfdG8gPSAicGFpciIsIHZhbHVlc190byA9ICJlbnRyb3B5X2RpZmZlcmVuY2UiKSAlPiUgCiAgZ3JvdXBfYnkocGFpcikgJT4lIAogIHN1bW1hcmlzZShlbnRyb3B5X2RpZmZlcmVuY2VfZGlzdHJpYnV0aW9uID0gbGlzdChlbnRyb3B5X2RpZmZlcmVuY2UpKSAlPiUgCiAgcmlnaHRfam9pbihwYWlyc19lbnRyb3B5X2RpZmYpICU+JSAKICBtdXRhdGUocF92YWx1ZSA9IG1hcDJfZGJsKGVudHJvcHlfZGlmZmVyZW5jZSwgZW50cm9weV9kaWZmZXJlbmNlX2Rpc3RyaWJ1dGlvbiwgfmVjZGYoLnkpKC54KSkpICU+JSAKICBtdXRhdGUocF9hZGogPSBwLmFkanVzdChwX3ZhbHVlLCBtZXRob2QgPSAiYm9uZmVycm9uaSIsIG4gPSBuKCkpKQplbnRyb3B5X2RpZmZfdGFibGUKYGBgCmBgYHtyLCBmaWcud2lkdGggPSA1LCBmaWcuaGVpZ2h0PTV9CiMgUGxvdCBkaXZlcnNpdHkgdmFsdWVzCmJyb29tOjp0aWR5KGRpdl9kaWFnKSAlPiUgCiAgbXV0YXRlKG5hbWVzID0gdG9sb3dlcihuYW1lcykpICU+JSAKICBnZ3Bsb3QoYWVzKHg9bmFtZXMsIHkgPSB4KSkrCiAgZ2VvbV9iYXIoc3RhdCA9ICJpZGVudGl0eSIpICsKICBnZ3B1YnI6OnN0YXRfcHZhbHVlX21hbnVhbCgKICAgIGVudHJvcHlfZGlmZl90YWJsZSAlPiUgc2VwYXJhdGUocGFpciwgaW50byA9IGMoImdyb3VwMSIsICJncm91cDIiKSwgc2VwID0gIlxcLiIpLAogICAgbGFiZWw9ICJwX2FkaiIsCiAgICB5LnBvc2l0aW9uID0gNywgc3RlcC5pbmNyZWFzZSA9IDAuMgogICkgKwogIHRoZW1lX2J3KCkgKwogIGxhYnMoeCA9ICIiLCB5ID0gIlNoYW5ub24gZGl2ZXJzaXR5IikgKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gNDUsIGhqdXN0ID0gMSkpCmBgYAoKCgojIyMgUGFpcndpc2UgcmFuZG9taXphdGlvbgoKYGBge3J9CnBhaXJzX2xpc3QgPC0gYXMubGlzdChhcy5kYXRhLmZyYW1lKGNvbWJuKHVuaXF1ZShkZiRjcml0ZXJpYSksMikpKQpuYW1lcyhwYWlyc19saXN0KSA8LSBzYXBwbHkocGFpcnNfbGlzdCwgcGFzdGUsIGNvbGxhcHNlID0gIi4iKQpwYWlyc19saXN0CmBgYAoKCgpgYGB7cn0KIyBDcmVhdGUgZnVuY3Rpb24gdGhhdCBzdWJzZXRzIGRmIG9mIGFsbCBkaWFnbm9zZXMgZm9yIGEgY3JpdGVyaWEgYmFzZWQgb24gYSBwYWlyIG9mIHR3byBjcml0ZXJpYSB0byBjb21wYXJlCiMgVGhlbiByYW5kb21pemVzIHRoZSBkaWFnbm9zZXMgYWNyb3NzIHRoZSB0d28gY3JpdGVyaWEsIAojIG1lYXN1cmVzIHRoZSByZXN1bHRpbmcgZGl2ZXJzaXR5IG9mIGVhY2ggY3JpdGVyaWEncyBkaWFnbm9zZXMKIyBhbmQgY2FsY3VsYXRlcyBhIGRpZmZlcmVuY2UgaW4gdGhlIGRpdmVyc2l0eSB2YWx1ZQojIHJhbmRvbWl6ZSBjYW4gYmUgc2V0IEYgdG8gb2J0YWluIG9yaWdpbmFsIGRpZmZlcmVuY2UgaW4gZGl2ZXJzaXR5IHcvbyByYW5kb21pemF0aW9uCiMgVG8gYmUgdXNlZCB3aXRoIHJlcGxpY2F0ZSB0byBnZW5lcmF0ZSBhIG51bGwgaHlwb3RoZXNpcyBkaXN0cmlidXRpb24gCnBlcm11dGF0aW9uX3BhaXJfaXRlcmF0aW9uIDwtIGZ1bmN0aW9uKHBhaXIsIHJhbmRvbWl6ZSA9IFQpewogIGRmIDwtIGRmICU+JSAKICAgIGZpbHRlcihjcml0ZXJpYSAlaW4lIHBhaXIpIAogIGlmIChyYW5kb21pemUgPT0gVCl7CiAgICBkZiA8LSBkZiAlPiUgbXV0YXRlKGRpYWdub3NpcyA9IHNhbXBsZShkaWFnbm9zaXMsIG4oKSkpCiAgfQogIGRpdiA8LSB2ZWdhbjo6ZGl2ZXJzaXR5KHRhYmxlKGRmJGNyaXRlcmlhLCBkZiRkaWFnbm9zaXMpKQogIHVubmFtZShkaXZbMV0gLSBkaXZbMl0pCn0KCnJlcGxpY2F0ZSgxMCxwZXJtdXRhdGlvbl9wYWlyX2l0ZXJhdGlvbihwYWlyc19saXN0W1sxXV0sIHJhbmRvbWl6ZSA9IFQpKQpyZXBsaWNhdGUoMTAscGVybXV0YXRpb25fcGFpcl9pdGVyYXRpb24ocGFpcnNfbGlzdFtbMV1dLCByYW5kb21pemUgPSBGKSkKYGBgCgpgYGB7cn0KIyAjIENyZWF0ZSBhIG51bGwgaHlwb3RoZXNpcyBkaXN0cmlidXRpb24gZm9yIHRoZSBkaWZmZXJlbmNlIGluIGRpdmVyc2l0eQojICMgRm9yIGFsbCBwYWlycyBvZiBjcml0ZXJpYQoKIyBzZXQuc2VlZCgxMjM0KQojIHBhaXJzX3Blcm11dGF0aW9uX2Rpc3RyaWJ1dGlvbnMgPC0gbGFwcGx5KHBhaXJzX2xpc3QsIGZ1bmN0aW9uKHgpIHJlcGxpY2F0ZSgxMDAwMCwgcGVybXV0YXRpb25fcGFpcl9pdGVyYXRpb24oeCkpKQojIG5hbWVzKHBhaXJzX3Blcm11dGF0aW9uX2Rpc3RyaWJ1dGlvbnMpIDwtIGxhcHBseShwYWlyc19saXN0LCBmdW5jdGlvbih4KSBwYXN0ZSh4LCBjb2xsYXBzZSA9ICItIikpCiMgZGF0YS5mcmFtZShwYWlyc19wZXJtdXRhdGlvbl9kaXN0cmlidXRpb25zKQojIHdyaXRlX2NzdihkYXRhLmZyYW1lKHBhaXJzX3Blcm11dGF0aW9uX2Rpc3RyaWJ1dGlvbnMpLCBoZXJlKCJkYXRhL2NyaXRlcmlhX3BhaXJ3aXNlX251bGxfZGlmZmVyZW5jZV9kaXN0cmlidXRpb25zLmNzdiIpKQoKcGFpcnNfcGVybXV0YXRpb25fZGlzdHJpYnV0aW9ucyA8LSByZWFkX2NzdihoZXJlKCJkYXRhL2NyaXRlcmlhX3BhaXJ3aXNlX251bGxfZGlmZmVyZW5jZV9kaXN0cmlidXRpb25zLmNzdiIpKQpwYWlyc19wZXJtdXRhdGlvbl9kaXN0cmlidXRpb25zCmBgYAoKCmBgYHtyfQojIENyZWF0ZSBhIGRmIG9mIHRoZSBhY3R1YWwgZGlmZmVyZW5jZSBpbiBkaXZlcnNpdHkgYmV0d2VlbiBwYWlycyBvZiBjcml0ZXJpYQpwYWlyc19lbnRyb3B5X2RpZmYgPC0gbGFwcGx5KHBhaXJzX2xpc3QsIHBlcm11dGF0aW9uX3BhaXJfaXRlcmF0aW9uLCByYW5kb21pemUgPSBGKQpuYW1lcyhwYWlyc19lbnRyb3B5X2RpZmYpIDwtIGxhcHBseShwYWlyc19saXN0LCBmdW5jdGlvbih4KSBwYXN0ZSh4LCBjb2xsYXBzZSA9ICItIikpCnBhaXJzX2VudHJvcHlfZGlmZiA8LSBkYXRhLmZyYW1lKHBhaXJzX2VudHJvcHlfZGlmZikgJT4lICAKICBwaXZvdF9sb25nZXIoZXZlcnl0aGluZygpLCBuYW1lc190byA9ICJwYWlyIiwgdmFsdWVzX3RvID0gImVudHJvcHlfZGlmZmVyZW5jZSIpICU+JSAKICBtdXRhdGUoZW50cm9weV9kaWZmZXJlbmNlID0gaWZlbHNlKGVudHJvcHlfZGlmZmVyZW5jZSA+PTAsIDAtZW50cm9weV9kaWZmZXJlbmNlLCBlbnRyb3B5X2RpZmZlcmVuY2UpKQpwYWlyc19lbnRyb3B5X2RpZmYKYGBgCmBgYHtyLCBmaWcud2lkdGg9MTIsIGZpZy5oZWlnaHQgPSA4fQojIENyZWF0ZSBwbG90cyBvZiBudWxsIGh5cG90aGVzaXMgZGlzdHJpYnV0aW9ucywgb3ZlcmxhaWQgd2l0aCBhY3R1YWwgZGl2ZXJzaXR5IGRpZmZlcmVuY2UKcmVhZF9jc3YoaGVyZSgiZGF0YS9hbGxfY3JpdGVyaWFfcGFpcndpc2VfbnVsbF9kaWZmZXJlbmNlX2Rpc3RyaWJ1dGlvbnMuY3N2IikpICU+JQogIGdyb3VwX2J5KHBhaXIpICU+JSAKICBuZXN0KCkgJT4lIAogIHNlcGFyYXRlKHBhaXIsIGludG8gPSBjKCJwMSIsInAyIiksIHNlcCA9ICJcXC4iLCByZW1vdmUgPSBGKSAlPiUKICByb3d3aXNlKCkgJT4lIG11dGF0ZShwYWlyPXBhc3RlKHNvcnQoYyhwMSxwMikpLCBjb2xsYXBzZSA9ICIuIikpICU+JSAKICBzZWxlY3QoIXBhaXI6cDIpICU+JSAKICB1bm5lc3QoZGF0YSkgJT4lIAogICMgcGl2b3RfbG9uZ2VyKGNvbHMgPSBldmVyeXRoaW5nKCksIG5hbWVzX3RvID0gInBhaXIiLCB2YWx1ZXNfdG8gPSAiZW50cm9weV9kaWZmZXJlbmNlIikgJT4lIAogIGdncGxvdChhZXMoeCA9IGVudHJvcHlfZGlmZmVyZW5jZSkpICsKICBnZW9tX2hpc3RvZ3JhbShiaW53aWR0aCA9IDAuMDEpKwogIGdlb21fdmxpbmUoZGF0YSA9IHBhaXJzX2VudHJvcHlfZGlmZiwgYWVzKHhpbnRlcmNlcHQgPSBlbnRyb3B5X2RpZmZlcmVuY2UpLCBjb2xvciA9ICJyZWQiKSsKICBmYWNldF93cmFwKH5wYWlyLCBzY2FsZXMgPSAiZnJlZSIpKwogIHRoZW1lX2NsYXNzaWMoKQpgYGAKCmBgYHtyfQojIENhbGN1bGF0ZSBwLXZhbHVlcyBmb3IgZGl2ZXJzaXR5IGRpZmZlcmVuY2VzIGJhc2VkIG9uIHRoZSBxdWFudGlsZSBvZiB0aGUgCiMgZGl2ZXJzaXR5IGRpZmZlcmVuY2UgdmFsdWUgY29tcGFyZWQgdG8gdGhlIG51bGwgaHlwb3RoZXNpcyBkaXZlcnNpdHkgZGlzdHJpYnV0aW9uIHZhbHVlcwojIFBhaXJ3aXNlIGNvbXBhcmlzb24gYWRqdXN0ZWQKZGF0YS5mcmFtZShwYWlyc19wZXJtdXRhdGlvbl9kaXN0cmlidXRpb25zKSAlPiUgCiAgcGl2b3RfbG9uZ2VyKGV2ZXJ5dGhpbmcoKSwgbmFtZXNfdG8gPSAicGFpciIsIHZhbHVlc190byA9ICJlbnRyb3B5X2RpZmZlcmVuY2UiKSAlPiUgCiAgZ3JvdXBfYnkocGFpcikgJT4lIAogIHN1bW1hcmlzZShlbnRyb3B5X2RpZmZlcmVuY2VfZGlzdHJpYnV0aW9uID0gbGlzdChlbnRyb3B5X2RpZmZlcmVuY2UpKSAlPiUgCiAgcmlnaHRfam9pbihwYWlyc19lbnRyb3B5X2RpZmYpICU+JSAKICBtdXRhdGUocF92YWx1ZSA9IG1hcDJfZGJsKGVudHJvcHlfZGlmZmVyZW5jZSwgZW50cm9weV9kaWZmZXJlbmNlX2Rpc3RyaWJ1dGlvbiwgfmVjZGYoLnkpKC54KSkpICU+JSAKICBtdXRhdGUocF9hZGogPSBwLmFkanVzdChwX3ZhbHVlLCBtZXRob2QgPSAiYm9uZmVycm9uaSIsIG4gPSBuKCkpKQpgYGAKCiMgUGxvdHMKCmBgYHtyfQoKYGBgCgo=